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Abstract 

We study a stochastic model proposed recently in the genetic literature to explain the 
heterogeneity of cell populations or of gene products. Cells are located in two colonies, whose 
sizes fluctuate as birth with migration processes in switching environment. We prove that 
there is a range of parameters where heterogeneity induces a larger mean fitness. 
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1 Introduction 



In P], the authors introduce a model for stochastic gene expression to study the heterogeneity 
of cell populations. They assume that the cells, or for example the product of some gene, can 
be in two distinct states or colonies. Let X(t) and Y(t) be the sizes of these colonies, which are 
here considered as birth with migration processes. We assume that the birth rates are either 71 
or 70 with A7 = 71 — 70 > 0, and that the associated migration rates k± and ko are such that 
ko ^ ki , that is cells located in the colony having the smaller birth rate migrate at a higher rate 
to the colony with the higher birth rate than the other way round. 

If the birth and migration rates are assigned once and for all to a corresponding colony (e.g. 
70 and ko to X, and 71 and k\ to Y), then the mean sizes no(t) = M(X(t)) and n\{t) = M(Y(t)) 
satisfy the pair of differential equations (see [Zj or [H] ) 

dno(t) , , . , . , ,,, 

dt = (70 - ko)no(t) + k 1 n 1 {t), 

= (-fi-h)ni(t) + k n (t). 

dt 

According to [OJ, we say that cells of the first colony represented by X(t) are unfit (they have 
the lower birth rates), and conversely that cells of the second colony represented by Y(t) are fit. 
The proportion of fit cells in the global population, y{t) = n\(t)/(n\{t) + no(t)), t ^ 0, satisfies 
the non-linear differential equation 

^ = k + (A 7 - k - h)y(t) - A 7 y(t) 2 . (2) 

Then, as t — > 00, y(t) — > yi, where \ ^ y\ ^ 1 follows directly from (j2j), see Section EJ 
This describes the equilibrium value of the proportion of fit cells in a non-changing environment. 
Fixing the values of the parameters ko, 70 and 71, we can ask for the value of ^ k± ^ ko which 
maximizes the proportion of fit cells, i.e. the equilibrium value of y{t): the optimal strategy is 
to keep all the fit cells in the fit state, that is to set their migration rate to zero, k\ = 0. This 
leads to y\ = 1, and thus the optimal solution would be a homogeneous population. 

Observations reveal however that most cell populations are not homogeneous; to explain 
this, the authors of 9 propose to introduce a small modification in the model by allowing 
environmental changes (for related questions in this context, see e.g. [5]), and show through 
Monte-Carlo simulations that the homogeneous solution k\ = is then not always optimal. The 
idea in their model is to allow the birth and migration rates to switch at random times from one 
colony to the other, so that cells in the fit colony become unfit and vice versa. If for example 
an environmental change occurs at some random time T\ > (To = 0), then the function fi(t) 
representing the proportion of fit cells solves © up to time Xi, and just after Ti, say at time 
T\ + 0, the fit cells corresponding to Y{t) become unfit and vice versa. The proportion of fit 
cells fi(t) is then switched to fi(T\ + 0) = 1 — f\(Ti). After T\, the random process {fi(t)}t^o 
solves (J2J) with initial data fi(T\ + 0) at time T\ + 0, until a new environmental change occurs, 
say at time T2 > T\. There is a new switch, and the process is again solution of @, until a new 
event occurs and so on. 

In [0], the fluctuations of the environment are modeled using a renewal process; the instants 
Ti, i ^ 0, are such that the sequence of random variables given by U = Ti — Tj_i, i ^ 1, is 

i.i.d. distributed according to some law /x on M + . The authors then use Monte-Carlo simulations 
to estimate the limiting value of the time averages along trajectories of the process fi(t), of the 
form 



1 [ ±N 

S N = 7^ A (*) da 
+ n Jo 



This limiting average value is denoted by Av(/i)fe 1 to express its dependency on the migration 
rate k\ < ko, when all the remaining parameters are fixed. Their simulations indicate that there 
is a range of parameters (ko not too large) such that 

Av(/i) fel>0 > Av(/i) fcl=0 , 

which means that heterogeneous populations are more adapted than homogeneous ones in a 
switching environment. 

In this paper, we study mathematically the limiting behavior of the stochastic process fi(t) 
and the associated time average Sn by giving its stationary measure, and we provide math- 
ematical formulas and numerical solutions, which might be of interest in practical laboratory 
experiments (see e.g. ). 

Our technique uses the process = /i(T^ — 0), Xq = /i(0), which is such that A^ +1 = 
Ptk+xO- — Afc), for some mapping ipt(x) (see Definition ^) . (Xk)k^o is a stochastic recursive 
Markov chain, and Sn can be expressed as an additive functional of the trajectory of (-XjOi^fc^iV- 
In Section 12 we recall a Theorem from Pj on the convergence of stochastic recursive chains, 
which applies in this setting. We give conditions ensuring the existence and uniqueness of a 
stationary measure tt, as well as geometric ergodicity. In Section |31 we consider the case where 
is exponential of parameter k > 0, and show that ir has a C°° density P with respect to 
Lebesgue measure. We furthermore prove in Theorem |2] that a multiple G of P solves a second 
order differential equation with weak singularities. Proposition ^ provides series expansions 
for P, which are necessary to derive properties of P near the singularities. In Section we 
show numerical solutions, using the series expansions of Proposition ^ to start the numerical 
integration. We provide an example where Av(/i)fc 1 >o > Av(/i)fc 1= o, which shows that it can 
be better to allow fit cells to migrate to the unfit state than to conserve all the fit cells in the 
fit state in such a switching environment. This is a regime where it is suitable for the colonies 
to anticipate bad hypothetical future events. 



2 Convergence of recursive chains 

We first give some basic results for the differential equation P|). The right hand side of 
can be factored into —Aj(y — yo)(y — y\), where yo = (A7 — ko — ki) — yd)/(2A r y) < 0, 
yi = ((A7- k - ki) + Vd)/(2Ai) > 0, and d = (A7 - k - h) 2 + Ak A7. Then k > h implies 
that 0<1 — ?/i<5< yx < 1, and that the derivative d/i(t)/dt is positive when fi(t) is in the 
interval [0, y\), negative in (yi,l], and it vanishes for fi(t) = y\. It is not hard to check that 
any realization of the trajectory {fi(t)}t^o, with initial data /i(0) G / = (1 — yi,yi) will remain 
forever in /, and that any trajectory starting in the interval I c = [0, 1] \ I will enter / after an 
almost surely finite time. (However, /i(0) = y\ implies fi(t) = y\.) We thus restrict our study 
to the interval /. 

Given t € M + , we define the mapping (ft : I — > I such that (ft(x) is the value of the solution 
of © at time t when starting at x £ I at time to = 0. Using separation of variables for we 
obtain the relation 

M x ) - VP x-yp 

T\ = exp(/3t), 3) 

y 1 -ift{x) y\-x 

where we set (3 = A r y(yi — yo). Given u G /, let 5t(u,y) denote the time interval the orbit of 
the dynamical system (j2J) needs to join u and y, y ^ u, when starting at time t = at u. Then 



(4) 



Definition 1 Given Xq = /i(0) G /, consider the Markov chain with values in I defined by 

where the sequence of random variables {tk}keN+ ^ s i-i-d. distributed according to some law \i 
on M + . This Markov chain describes the evolution of fi(Tk — 0), at the instants just before the 
switches, with T^ +1 — T& = . 

We first recall and adapt results of [3] on the convergence of such Markov chains, also called 
stochastic recursive chains, see e.g. pQ. The general setting is described by a complete separable 
metric space (S,p), the set of values taken by the Markov chain, a family of mappings fg : 
S — * S, indexed by parameters 9 living in some parameter space 0, and a probability measure 
p on Q. Given an i.i.d. sequence of random elements 9 n , n ^ 1, of law p, we can consider the 
Markov chain (X n ) n£ ^ given by X n+ \ = fe n+1 (X n ). The following Theorem gives conditions 
for the existence and uniqueness of a stationary measure (Theorem 1.1 of [2]) • In what follows, 
p( n \x, dy) denotes the law of the Markov chain X n and /o[-P( n )(x, •), vr] is the Prokhorov metric, 
see below. 

Theorem 1 Assume that the family of functions fg, 9 G G is Lipschitz with 

p{fe{x),fe{y))^K e p{x,y), x,y G S, 
V# G G. Assume furthermore that 

J K 6 fj,(d9) < oo, J p(fg(x ),x )p(d9) < oo, (5) 

for some xq G S, and that 

f \n{K e )p{d9) < 0. (6) 

Then 



• The Markov chain has a unique stationary distribution tt, 

• p[p( n \x, •), tt] ^ A x r n , for constants A x and r with < A x < oo and < r < 1; this 
bound holds for all times n and all starting positions x, 

• the constant r does not depend on n or x; the constant A x does not depend on n, and 
A x < a + bp(x, xq), where < a, b < oo. 

In our setting, S is given by / and the parameter set Q is just R + . The Prokhorov distance 
d n := p{p( n ^ (Xq, •), vr] is the infimum of the 5 > such that 

P^{Xq,C) <7r(Cs)+5 and tt(C) < P (n) (*o, C s ) + 6, (7) 

where C runs over the Borel sets of I and, for given C G B(J), Cs denotes the set of points of 
/ whose distance from C is less than 5 (see Section 5.1 of (Sj). Condition (jEJ) means that the 
functions fg are contractions in the average. We first express this condition in our setting: for 
t G Q = M + and u £ I = S, the mapping ipt{u) is given explicitly by 

, x _ 2/o(2/i ~u) +yi(n-y )exp(/3t) 
MU) ~ yi -u + (u-yQ)eMPt) ■ (8) 

Setting ft(x) = (pt(l — x), we obtain 



Lemma 1 For all t 6 



d , , x (yi - yo) 2 exp(pt) 

:ft{x) - 



dx (yi — 1 + x + (1 — x — y ) exp(/?i)) 2 ' 

^ ,d, ni (2/1 - y ) 2 exp(^) 

^ :=SUp — / t (x) = — — — r T^TVTa- 

xeI dx (2yi - 1 + (1 - yx - y ) exp(/3t)) 2 

// /i is exponential of parameter k > 0, and a = i/ien i/ie conditions given in |2J) are 

satisfied, and 

f K exp(-Kt)]n(K t )dt = -a-2z [°° eXp( ~ (1 ± dt, 
Jr+ Jo 1 - z exp(-i) 

where we set z = —(2yi — 1)/(1 — yi — yo) < 0. Condition |^ is t/ius satisfied if 

f°° exp(-(l + a)t) , 

- a - 2z / PV v ^Vdt < 0. 9) 

7 l-zexp(-i) 

Remark 1 When \z\ ^ 1, the integral j ° o (exp(— (1 + a)t))/(l — zexp(— t))dt is the Lerch Phi 
function <&(z, s, v) = Yln>o( v + n )~ Szn , with s = 1 and i> = 1 + a, and is also equal to Gauss's 
Hypergeometric function 2Fi(l,l + a; 2 + a; z)/(l + a) (see e.g. chap. 1.11). 

Proof: Taking the derivative of Q with respect to u, we obtain 

d , . (yi- 2/o) 2 exp(/3t) 

—(pt(u) = 2 

du ( yi _ u + ( u - yo ) exp(pt)) 

and thus 

-// \ _ __d_ r-, _ \- (yi - yo) 2 exp(ffl) 

Mj " dn^ U X) ~ ( yi -l + x + (l- x - yo )eMPt)) 2 
as required. The expression for K t follows from direct computation. 



3 Convergence to stationarity in Poissonian environments 

Assume that /i is exponential of parameter k > 0. We will see in the sequel that the stationary 
measure ir has, under some conditions, a density P(y) such that with Q(y) = ((y — Vo)/(yi — y)) a , 
where a = k/(3, the function G(y) = P(y)Q(y)(y — yo)(yi — y) satisfies the differential equation 

G"(y) + U(y)G'(y) + V(y)G(y) = 0, (10) 

where y = 1 - yo, yi = 1 - yi, 

£/(y) = — - — + — — , (11) 

y - yo y-yi y-yi y - yo 



and 



y (y ) = <* 2 (2/i-2/o) 2 (12) 

(y - yo)(y - yi)(y - yo)(y - m) 



The following proposition will therefore be useful: 

Proposition 1 The solutions of the second order homogeneous linear differential equation U(J\) 
are analytic on the interval I = (yi,yi). Two fundamental solutions G\(y), G^iy) are 



G\{y) = (y — yi) a W\{y) , where W\{y) is analytic on (yi — 5,yi) for some 5 > and with 

wiG/i) = i. 



G 2 (y) 



W 2 (y), if a £ Z, 

W 2 (y) + CG 1 (y)Hy-y l ), if a G Z, 
wif/j Vt^Cy) analytic on (jji — 6,yi) for some 5 > 0, ^2(2/1) = 1 and C G 



Another set of two fundamental solutions G\{y), G 2 (y) is 

• G\{y) = (y\ — y) l ~ a W\{y), where W\{y) is analytic on (yi,yi + 5) for some 5 > and 
with Wi(y\) = 1. 



W 2 (y), ifagZ, 
W 2 (y) + CGi(y)ln( yi -y), ifa€Z, 
u>i£/i ^2(2/) analytic on (yi,yi + 5) /or some <5 > 0, W 2 (y\) = 1 and C G 



G 2 (y) 



In the appendix, we show this result for completeness, and also how these fundamental solutions 
can be computed by series expansion about y± and y\ respectively. 

Theorem 2 Assume that 

f°° exp(-(l + a)t) , 
-a - 2z / y y — r-^r-dt < 0, 
J l-zexp(-i) 

where z = — (2yi — 1)/(1 — yi — yo) < 0. Then the Markov chain X^ from Definition^ with 
initial data Xq G / = (1 — yi,y%) has a unique stationary distribution tt of C°° density 

p , s = Q{y)-\y - yi) a w l (y)/( yi - y )/( y - yo ) 

[V> fj Q{z)~\z - m) a W 1 (z)/(y 1 - y)/(y - y ) dz 

Here, Q{y) = (jj^fj > where a = k/(3, W\(y) is the analytic function on {y\ — 5,y\) with 

W(yi) = 1, such that G%(y) = (y — yi) a W\{y) is a solution of the differential equation H0\) . In 
the neighborhood of y = y\, this solution is such that < lim^^ W\{y) < +00. Finally, the 
behavior of the density P near y\ is given by (y\ — y) a ~ l , and thus converges when a ^ 1 and 
diverges toward +00 when a < 1. Let f(x) = x and g{x) = ln((x — 1 + yi)/(yi — x)) be defined 
on I. Then g G 1/(7, B(J), tt) with 

= Vo + 7^(0). (13) 

Remark 2 Relation ilty) will be useful when considering time averages for Monte- Carlo simu- 
lations, see Section^ 

Proof: The existence and uniqueness of the stationary measure follows from Theorem ^ and 
Lemma m Let Y be a random variable of law tt, and let T be exponential of parameter k > 0, 
independent of Y. In the stationary regime, Y =c ^r(l — Y). Let F(y) = P(Y < y). Then 



F (y) = / n(dv)Kexp(-Kt)I(ip t (l - v) < y)dt, 

JlxM+ 

where !(•) denotes the indicator function. For given y G /, the time variable t is restricted to 
the interval ^ t < 6t(yx, y) , see (@J. Thus 

fSt(l-yi,y) /• 

F{y) = / Kexp(-Kt) / 7r(du)% t (l - v) < y)dt. 



For given t in this interval, the set of v E I with tpt(l — v) < y is given by 



r ^ r i ^ 2/1(2/ - Vo) + exp(/3i)(yi - y)y 

{v £ i; 1 — v < r }. 

y-yo + exp(/ft)(yi -y) 

It follows that 7r(du)I(^(l— v) < y) = 1 — F(l — u), where we set u = (2/1(3/ — 2/o)+exp(/3i)(yi — 
y)Vo)/(y — Vo + exp(/3i)(yx — y)), with t = 5i(it, y). We make the change of variable i = <5i(u, y) 
with 

d* yx - yp 

du P(yi-u)(u-y )' 

Then 

KI Vy-yo; } x _ yx ( yi -u)( u - yo ) \ yi -J V V " 

This is a fixed point equation for the distribution function F. We use it for proving that the 
probability measure tt has a C°° density on the interval I. First notice that F is monotonically 
increasing and integrable on /. The above relation then shows that F is continuous on /. Using 
again this argument recursively, one sees that F is the integral of a continuous function and 
is therefore differentiable, with a continuous derivative. The C°° differentiability is obtained 
by iterating this argument. Let P be the C°° density of ir with respect to Lebesgue measure. 
Our strategy runs as follows: We use the fixed point relation to show that a multiple G of P 
satisfies a second order differential equation, which has only weak singularities, and then deduce 
properties of P with the help of Proposition^ 

For given v £ I, the time variable t is restricted to the interval 

< t ^ 8t(u, y) = ln((y - y )(yi - y)/(yi -y){u- yo))/P, 

where u = 1 — v (see HJ . It follows that 

ryi fSt(u,y) 

>l-y 



ryi fSt(u,y) 
F(y) = / P(v)dv / Kexp(-Kt)dt, 
Ji-v Jo 



with 



"(y) = — l = / -Puudu/cexpf — Kot(u, y)j 

dy Ji-y dy 



ryi 

a / dvP(v) 

Jl-V 



Q(u) (yi - yo) 



i- y ' Q(y) (y- 2/0X2/1 - vY 

where we set Q(y) = ((y - yo)/(yi - y)) Q - Using u = 1 — v and setting G(y) = P{y)Q{y){y - 
2/o) (2/1 - y), one gets 

ry 

G(y)= G(l-u)R(u)H(u)du, (14) 
Ji-vi 

where R(u) = aQ(u)Q(l — u) _1 is such that R(l — u) = a 2 /R(u), and H(u) = (yi — yo)/(yi — 
1 + it)/(l — u — yo)- Taking the derivative gives 

G'(y) = G(l- y )R(y)H(y), (15) 

or 

G(l - y) = G'(y) J R(y)- 1 i/(y)- 1 = a- 2 G'(y)i?(l - y)/H(y). 
Taking a second derivative then gives 

G"(y) + A ln(fc^)G'(y) + « 2 F(y)^(l - y)G(y) = 0. 



and simplifying the terms leads to (jlOj) . We see that R(u)H{u) ~ (u—l + as u — > \ — y\. 

The exponents associated with the fundamental solutions are /O = or a in the neighborhood of 
y = 1 — yi and p' = or 1 - a near y = y\. 

Assume first that a N + . We first check the behavior of G in a neighborhood of y = y\. Set 
y = y\ +e, e > 0, with 1 — y = y\ — e. Proposition^shows that G is a linear combination G(y) = 
Ae a W 1 (y)+BW 2 (y), for constants A, B £ 1 Similarly, G(l-y) = Ae 1 - a W l {l-y)+BW 2 {l-y), 
for real constants A and .B. As e — ► 0, the right hand side of l|14|l behaves like e a G{y\ — e) — » 0. 
Suppose that B / 0. Then G(y) ~ 7^ 0, and (fUjl can't be satisfied. One must thus 

have B = 0, so that G(y) = Ae a W 1 (y). When a > 1, (EI) implies that A = 0. It follows that, 
for arbitrary a > 0, lim^^ G{y) = BW2(yi) / 0, and that G(yi + e) ~ Ae a Wi(yi), e — > 0, as 
required. The corresponding result for P follows. 

Suppose that a £ N + . The right hand side of (|14|) behaves like 

F(e) := e a (Ae 1 - a W 1 (y 1 ) + B(W 2 (yi) + Ce^WM ln(e))), 

with F(e) — ► as e — ► 0, and G(yi + e) behaves like 

F(e) := Ae a W x [y{) + B(W 2 (yi) + Ce»W x {vx) M?))- 

One has F(e) ~ W 2 (yi), £ ^ 0, when S ^ and F(e) ~ ie a #i(yi), when 5 = 0. ED 
shows that necessarily B = 0. Suppose that a = 1. Then one must have BC = 0, implying 
the existence of the limit lim^^ G(y) ^ 0. When a > 1, A = 0, B ^ 0, and lim^^ G(y) = 
^^2(2/1)1 as required. 

Finally, we check the identity (|13j) . First 5 G L 1 (/, B(J), 7r) follows from the behavior of the 
density P at the boundaries of /, as described above. Next, 

E^g) = [ ln( y ~ 1 + yi )P(y)dy, 

J 1 yi-y 

where J := j I ln(y — 1 + y\)P(y)dy is such that 

J = / Hy-l+ yi )G(y)Q(y)- l H{1 ~ y) dy 

J 1 yi - yo 

J ln(yi - «)G(1 - u)Q(l - u) _1 i?(u)d 
-G(l - u)R(u)H(u)du 
G'(u)du 



u 



yi - yo 

1 f ln(yi - u) 



"(yi - yo) J 1 Q(u) 

1 /" ln(yi - u) 



«(yi - yo) J 1 Q(u) 



"(yi - yo) v Q(u) 7/ v QO) 

1 /" G(u) (tt - y 



' ^(n)i%^ yi - /G(n)fi% T ^)dn 



"(yi - yo 

where we use (EH). It follows that 



f G(u) (w- yo) , . /\ / u 
J 1 Q{u) (yi - u){u - y ) 7/ 



iMy) = 7 — ~ — r E -(/) 



«(yi - 2/0) <*(yi - y ) ' 

proving (|T3|) since a = K,/(Aj(yi — yo))- 

Corollary 1 Assume that condition holds. Then, as t —* +00, £/ie law 0/ i/ie stochastic 
process fi(t), t ^ 0, /i(0) G J, converges toward the stationary measure n of density P of the 
Markov Chain X^. 



Proof: Given t € M + , let t* be the last renewal time before t, and set S* = t — t*. When 
the length of the overlapping random interval is exponential, S* is also exponential. In the 
stationary regime, or equivalently for large t, one has the identity in law fi(t) =c <PS*0- ~ X), 
where X is distributed according to it, and the result follows. 



4 Time averages 

When the conclusions of Theorem ^ hold, the chain X^ has a unique stationary probability 
measure it, and Y^k=i 9(Xk)/n converges a.s. toward the expectation of g under it, for any 
function g in L l {I, B(/),7r), (see e.g. 0). In [5], the authors use Monte-Carlo methods based 
on the process fi(t), t ^ 0, to estimate the mean fitness by considering the time average 



1 



S N = ^J o /lOOda, (16) 
where N is a fixed number of renewal periods. 

Lemma 2 Lei ./V G N + . Given a realization = To < T\ < • • • < of the renewal process, we 
have 

i [ TN f( , A , (yi-yo) , ^ ^-i-q-j/i) ^ 
^y /l(s)ds=yo+ ^r^ ln (n yi _ Xi )• ( 17 ) 



Proof: Consider the integrals 



/i(s)ds, 



i-l 



where /i(Tj__i + 0) = 1 - and /i(T, - 0) = X { . The value of y(s) = f\(Ji-\ + s), 

s E (0, Tj — Tj_i) is given implicitly by Q; Therefore 

/ \ 2/o(yi - u) + yx{u - yo ) exp(ps) 
yx-u+[u- y Q ) exp(ps) 

where we set u = 1 — Xi~\, and thus, after a longer but not difficult computation, one obtains 

f Tl t , s A frr r \ i Ux Uo 1 fyi ~ u + (u - y ) exp(/3(Tj - Tj_i)) 
/ /i(s)ds = y (Ti -Ti_i) H — In 

and the result follows, since 

y x - u + (u - yo) exp(/3(Ti - T H )) = (yi-it)(l + - — — exp(/3^)) 

yi - it 

(yi - (1 - X t ^))( yi - y ) 

yi-Xi 

Theorem 3 Suppose that /j, is exponential of parameter k > 0, and assume fQ). Let f(x) = x 
and g(x) = ln((x — 1 + yi)/(yi — x)) be defined on I. Then 

1 [ Tn k 
lim—/ fi(s)ds = y Q + —Kn(g)=Rn(f), a.s. 
N-^oo In Jo A 7 



Proof: From equation (|17|). we obtain 



1 [ N * r \i , ivi-yo) , f X -l + yi (yx-yo 



his) ds = y + ln( " u V ) + > 



Ijv Jo " /^iv yi - ^Ov /Wj 



iY 



=1 



As Tjv is a renewal process with exponential inter arrival times of parameter k, it follows that 
Tn/N converges a.s. toward 1/k. Next, g £ B(J), 7r) follows from the behavior of the 

density P at the boundaries of J, as described in Theorem [21 From Proposition ^ and Theorem 
[21 the behavior of P in the neighborhood of y = 1 — y\ is given by (y — 1 + yi) pi where p\ = a 
and by (y\ — yY 2+a ^ 1 in the neighborhood of y = yi, where /?2 = 0. The Markov chain is 
geometrically ergodic, and thus the last term converges a.s. toward (K/(A7))E 7r (g). We finally 
check that ln(yi — Xn)/N converges a.s. toward 0. Given e > 0, consider the probability 

P(| ln(y! -X N )\> Ne) = P(ln( yi - X N ) < -Ne) 

= P(X N > yi - exp(-iVe)) = pW( Xo ,A N ), 

where An = {x > y\ — exp(— Ne)}. Using the behavior of P in the neighborhood of y = 
yi, one gets that tt(An) ^ M(exp(— eN)) P2+a , for some positive constant M. Let '■= 
|P (iV) (X , A N ) - ir(A N )\, and let d N be the Prokhorov distance defined in 0. If ir(A N ) ^ 
PW(X ,^), then 1N ^ tt{A n ). If tt^jv) ^ P^ N \X ,A N ), one has P^(I 0) 4) < 
ir((AN)d N ) + dN, and it follows that 

1 N = pW(X ,A N )-n{A N ) < 7r((A N ) dN ) - tt(A n ) + d N 

rVi —cxp(-eN) 

= / P(y)dy + d N 

^ d N + P>(exp(-eA0 P2+Q - (d N + exp(-eiV)) p2+Q ), 

for some positive constant D > 0. Theorem ^ gives that 

P(\H yi -X N )\>eN) ^ \P^\X ,A N )-n(A N )\+7r(A N )^h(X )X N , 

for some bounded function h and a positive number < A < 1. The result then follows from 
the Borel-Cantelli Lemma. The last identity is of Theorem [21 



5 Numerical Examples 

We now compute the density P given in Theorem[2]numericalry. To do so, we solve the differential 
equation Q1U|) numerically, starting in the neighborhood of the singular point y = y\ = 1 — y%. 
Proposition ^ an d Theorem [21 show that lim^-^ P(y) = 0, and that the first derivative of P 
behaves like (y — yi) a_1 , which goes to +oo when a < 1. We start the numerical solution at the 
point y = yi + e, where e > is small, and use the initial conditions G{y\ + e) and G"(yi + e) 
from the series expansions given in Proposition^ In addition, we use the numerical integration 
procedure to compute the integral to scale the density P, by adding an additional ordinary 
differential equation to (jlUj) . 

We show in Figures ^ to [3] the results obtained for five different sets of parameters. In all 
the figures, we show the computed solution G of the differential equation (jTU|) in dashed, the 
computed density P as a solid line, and the results of a Monte-Carlo simulation with lOO'OOO 
samples as circles. The density from the theory and the Monte-Carlo simulations agree very 
well. It is interesting to see in Figures ^ and [21 the variety of densities that can be generated 
by this simple model. Figure [21 contains a case where increasing k\ increases the overall fitness 
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Figure 1: Density P on the left when k = 10, A7 = 1, ko = 0.4, k\ = 0.05, and on the right 
when k = 1.5, A7 = 1, k = 0.1, ki = 0.05. 
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Figure 2: Density P on the left when k = 10, A7 = 9, ko = 0.1, k\ = 0. Av(fi)k 1= o = 
0.553274111, and on the right when k = 10, A7 = 9, k = 0.1, k x = 0.05. Av(f{) kl=om = 
0.55672212. Clearly the average fitness is larger when k\ = 0.05 than when k\ = 0. 
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Figure 3: Density P when k = 5, A7 = 3, &o = 3, k\ = 1, a case where a < 1. 

of the population. Figure El finally shows a case where a < 1. We note that the numerical 
integration out of the singularity can be challenging. In particular, for the first case in Figure 
^ the standard ode45 from Matlab needed very small absolute tolerances to succeed with the 
integration for e < le — 2. A more robust method turned out to be DOPRI853, see [Sj. 
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6 Appendix 

In this appendix we show for completeness the proof of Proposition ^ an d describe a method 
how to solve the differential equation (|TT)1) (see also pp. 317-321). This equation is of the 
form 

G"(y) + U(y)G'(y) + V(y)G(y) = 0, 

where the functions U{y) and V{y) are meromorphic in the complex plane with four poles of 
order one at yo < Vi '■= 1 — Hi < Vi < Vo '■= 1 — 2/o- The solutions are therefore analytic in 
the open disc of radius (yi — y*i)/2 centered at 1/2. We look for real solutions in the interval 
I = (2/1,2/1). In order to simplify calculations, we use the variable transformation 

y = vi + {yi - yi)z, z = - — -z- (18) 

yi - yi 

and set g(z) := G(y). With this transformation, the differential equation (|10j) becomes 

g"(z)+u(z)g'{z) + v(z)g(z) = 0, (19) 
where u and v have four poles of order one at the points —b <0<l<l + 6 with b = 

(m - yo)/(yi - m)- 

, . 1 — a a a a + 1 , s ail + b) 2 

u{z) = — — + — -r-— -r + - t^—tt , v(z)- 



z-1 z + b z -(l + by w z(l - z){z + b)(l + b - z) 
We can therefore rewrite this equation as 

g"{z) + ^-g'{z) + ^g(z) = 0, (20) 
where h(z) and k(z) are analytic in the disc of radius min{l, b} centered at 0: 

oo oo 

h(z) = Y,a n Z n , k{z)=Y J Pn* n - 
n=0 n=0 

Multiplying the equation ()2U|) by z 2 we get an equivalent equation which can be written as 

L(g) := (/x 2 2 D 2 + MhMjg D + fi k )(g) = 0, (21) 

where D denotes differentiation and fj,f multiplication by a function f(z). Looking for solutions 
of the form 

oo 

g(z) = z p w(z), w(z) = 1 + ^ w n z n , 

n=l 

we may identify the function g{z) with the infinite row [w] = [1, wi, W2, W3, ■ ■ ■] and write ()21|) 
in matrix form: 

[L/>]M T = 0. (22) 
If we write L as L = (fJ, z D + /i/ l _ 1 )/i 2 D + /i^, we get the lower triangular matrix 



[1/ 



p(p + a - 1) + Po 

pai+Pi (p + l)(p + a ) +Po 

pa 2 + P 2 (p+l)oi+/3i (p + 2)(p + l + a ) + O 



A solution [w] = [1, w%, W2, ■ ■ ■] of the linear system (|22|) exists if and only if Lq = 0. This is the 
so-called indicial equation for p. From now on we shall no longer treat the general case but only 
the case corresponding to our differential equation (|T9|) . In this case ao = I — a and /3o = 0. 
So the indicial equation is p(p — a) = and yields the two characteristic exponents p\ = a and 
P2 = 0. We shall write instead of Lfj. 

For p = pi, the solution [it/ 1 )] = [1, w±, ■> ■ ■ ■] ma y be calculated by the recursion scheme 

4 X) = 1, = TT~ IE Lij^A for n > 1. 

nn \j=o J 

With these coefficients Wn \ the function 



) z n 



n=l 



is a solution of ()19|) . From the general theory of linear differential equations in the complex 
plane it follows that g\ is analytic in the disc of radius 1/2 centered at 1/2, but the power series 
for w\(z) might have a convergence radius < S < 1. 

If a is not an integer, another solution g 2 (z), linearly independent of gx(z), can be obtained 
in the same way from p = p2 = 0. If, however, a is an integer, the corresponding matrix 
has the entry L^ n = for n = a, and we look in this case for a solution 52 i z ) of the form 
92(2) = 1 + En^i^n z " + Cgi(z) Inz. As 51 is a solution, the terms in L(g 2 ) containing Inz 
cancel and the function w^{z) = 1 + J2 n >i ^n' 2 " must satisfy the equation 



L( w ( 2 )) = -C(2p z B + p h -i){gi). 

Identifying w^ 2 \z) with the infinite row [w^\ = [1, . . .], we can write this in matrix 

form 

[L 2 ][ W ^ = -C[v 1 ,v 2 ,...f. (23) 

For the right-hand side one checks easily that Vj = for j = 0, . . . , a — 1 and v a = a. Therefore 
we can resolve the inhomogeneous linear system (|23|) in the following way: 

1. We determine w^p for j ^ a in the same way as Wj . 

2. We set w a ■= and determine the constant C by the equation Y^j=o L a jWj(2) = —Cv a . 

(2) 

3. We determine the coefficients Wn for n > a by the recursion formula 



nn \ j=0 J 



W M = — I Cv n . + I for n > a + 1. 



We shall not go into further details, for example present concrete formulas expressing the v n 
by the Wn \ because we don't really need the solution g 2 of (|2H|) in our case, as we have shown 
in the proof of Theorem [2 

Using the variable transformation (|18|) we get the solutions Gj(y) of the original differential 
equation pi)|). in particular 

~ \ / 00 (i) \ 

y-yi\ / ~ \ d TT7" ^ / - \a / -, , ^ \ 



Gi(y) = (yi - (^-4^) = (y - yi) a Wi(v) = (y - yi) Q ( 1 + £ 



^(yi-yi)" y 



In order to find fundamental solutions near the singularity y\ , we can apply the same method 
once more, but using the variable transformation 



y = yi-{yi-yi)z, z = — . 

y\ - yi 

One easily checks that in this case the indicial equation is p(p + a — 1) = and that therefore 
the two characteristic exponents at y\ are p'-y = 1 — a and p' 2 = 0. We obtain thus the second 
fundamental system of solutions G\{y) and G2(y). 
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